{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/Qiskit/qiskit-tutorials/master/images/qiskit-heading.png\" alt=\"Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook\" width=\"500 px\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _*Hamiltonian Parameters*_ \n",
    "\n",
    "\n",
    "***\n",
    "### Contributors\n",
    "David McKay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The basic architecture of a quantum device for superconducting transmon qubits [[1](#ref1)] is that of several qubits coupled together by resonant bus structures. This can be described by the Hamiltonian \n",
    "$$H = \\sum_{i=0}^{N_{\\text{qubits}}} \\left(\\omega_i \\hat{n}_{i} + \\alpha_i \\frac{\\hat{n}_i}{2} (\\hat{n}_i-1) \\right) + \\sum_{j=0}^{N_{\\text{bus}}} \\omega_{\\text{bus},i} \\hat{n}_{\\text{bus},j}+ \\sum_{i,j} g_{ij} \\left(\\hat{a}_i \\hat{b}_j^{\\dagger} + \\hat{a}_i^{\\dagger} \\hat{b}_j\\right) $$\n",
    "\n",
    "where $\\omega_i$ is the $0\\rightarrow 1$ frequency of qubit $i$, $\\alpha_i$ is the anharmonicity of the transmon qubit $i$, $\\omega_{\\text{bus},j}$ is the frequency of bus $j$ and $g_{ij}$ is the coupling of qubit $i$ to bus $j$. In a planar architure there will only be a few non-zero values of $g_{ij}$. \n",
    "\n",
    "**Contents**\n",
    "\n",
    "[Dressed Basis](#sect1)\n",
    "\n",
    "[ZZ](#sect2)\n",
    "\n",
    "[Dispersive Shifts](#sect3)\n",
    "\n",
    "\n",
    "### References\n",
    "\n",
    "[1]<a id=\"ref1\"></a> Jens Koch, Terri M. Yu, Jay Gambetta, A. A. Houck, D. I. Schuster, J. Majer, Alexandre Blais, M. H. Devoret, S. M. Girvin and R. J. Schoelkopf. Charge insensitive qubit design derived from the Cooper pair box. https://arxiv.org/abs/cond-mat/0703002\n",
    "\n",
    "[2]<a id=\"ref2\"></a> Jerry Chow. Quantum Information Processing with Superconducting Qubits. https://rsl.yale.edu/sites/default/files/files/RSL_Theses/jmcthesis.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code imports and basis functions\n",
    "=============="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-18T15:13:54.475305Z",
     "start_time": "2018-12-18T15:13:53.026353Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "# qubit frequency (fq) dressed by resonance fres\n",
    "# interacting with g\n",
    "def qubit_dressed(fq, fres, g):\n",
    "    return fq + g**2/(fq-fres)\n",
    "\n",
    "# exchange interaction between qubits coupled to a \n",
    "# common bus\n",
    "def J(fq1, fq2, fbus, g1, g2):\n",
    "    d1 = (fq1-fbus)\n",
    "    d2 = (fq2-fbus)\n",
    "    return g1*g2*(d1+d2)/2/d1/d2\n",
    "\n",
    "# state dependent cavity shift\n",
    "def Chi(fq1, fbus, g, alpha):\n",
    "    d1 = (fq1-fbus)\n",
    "    d2 = (fq1+alpha-fbus)\n",
    "    return g**2*alpha/d1/d2\n",
    "\n",
    "# build an operator from a list\n",
    "# using kronecker product\n",
    "def multiqop(curop=None, oplist=None):\n",
    "    if curop is None or len(curop)==0:\n",
    "        return multiqop(oplist[0],oplist[1:])\n",
    "    if oplist is None or len(oplist)==0:\n",
    "        return curop\n",
    "    \n",
    "    return multiqop(np.kron(curop,oplist[0]),oplist[1:])\n",
    "\n",
    "def nop(n=1):\n",
    "    return np.diag(np.arange(n+1,dtype=complex))\n",
    "\n",
    "def aop(n=1):\n",
    "    aop_tmp = np.zeros([n+1,n+1],dtype=complex)\n",
    "    for i in range(n):\n",
    "        aop_tmp[i,i+1] = np.sqrt(i+1)\n",
    "    return aop_tmp\n",
    "\n",
    "def iop(n=1):\n",
    "    return np.eye(n+1, dtype=complex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sect1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dressed Basis\n",
    "\n",
    "The first assumption made in these systems is that the bus frequency and qubit frequencies are far apart and that the bus is not thermally occupied. Under these assumptions we can rewrite the Hamiltonian in a basis with the bus degrees of freedom removed and an effective coupling between qubits (rewriting without the anharmoncity for simplicity),\n",
    "$$H = \\sum_{i=0}^{N_{\\text{qubits}}} \\gamma_i \\hat{n}_{i} +  \\sum_{i,j<i} J_{ij} \\left(\\hat{a}_i \\hat{a}_j^{\\dagger} + \\hat{a}_i^{\\dagger} \\hat{a}_j\\right) $$\n",
    "\n",
    "where\n",
    "$$\\gamma_i = \\left[\\omega_i+\\sum_{j}\\frac{g_{ij}^2}{\\omega_i-\\omega_{\\text{bus},j}}\\right] $$\n",
    "\n",
    "Noting that the operator $\\hat{n}_i, \\hat{a}_i$ changed as well. Assuming a pair of qubits are coupled through only a single bus, the exchange term is expressed as ,\n",
    "$$ J_{ij}=\\frac{g_{ik} g_{jk} (\\omega_i+\\omega_j-2\\omega_{\\text{bus},k})}{2(\\omega_i-\\omega_{\\text{bus},k})(\\omega_j-\\omega_{\\text{bus},k})} $$\n",
    "\n",
    "We can further diagonalize,\n",
    "\n",
    "$$H = \\sum_{i=0}^{N_{\\text{qubits}}} \\lambda_i \\hat{\\tilde{n}}_{i} $$\n",
    "\n",
    "where \n",
    "\n",
    "$$\\lambda_i = \\gamma_i + \\frac{J^2}{\\gamma_i-\\gamma_j}$$ \n",
    "\n",
    "and so eventual qubits (the two-level systems in the final diagonalized Hamiltonian) are not the same as the bare qubits in our \"lab\" Hamiltonian. \n",
    "\n",
    "**Code**\n",
    "\n",
    "The code below numerically constructs the Hamiltonian for 2 qubits (infinite anharmonicity) coupled to a single bus. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Hamiltonian with n qubits coupled to a single bus\n",
    "def Hq(fq, fbus, g, alpha=None, qlevels=1, buslevels=1):\n",
    "    \n",
    "    oplist_id = [iop(qlevels) for e in range(len(fq))]\n",
    "    oplist_id.append(iop(buslevels))\n",
    "    \n",
    "    if alpha is None:\n",
    "        alpha = [0 for e in range(len(fq))]\n",
    "    \n",
    "    H = 0*multiqop(oplist=oplist_id)\n",
    "    \n",
    "    for qind, fqi in enumerate(fq):\n",
    "        oplist_i = oplist_id.copy()\n",
    "        oplist_i[qind] = nop(qlevels) \n",
    "        \n",
    "        #add qubit frequency\n",
    "        nqi_fullop = multiqop(oplist=oplist_i)\n",
    "        H += fqi*nqi_fullop\n",
    "        \n",
    "        #add anharmonicity\n",
    "        H += 0.5*alpha[qind]*nqi_fullop*(nqi_fullop-1)\n",
    "        \n",
    "        oplist_g0 = oplist_id.copy()\n",
    "        oplist_g1 = oplist_id.copy()\n",
    "        oplist_g0[qind] = aop(qlevels) \n",
    "        oplist_g1[qind] = aop(qlevels).transpose() \n",
    "        oplist_g0[-1] = aop(buslevels).transpose() \n",
    "        oplist_g1[-1] = aop(buslevels)\n",
    "        \n",
    "        #add bus interaction\n",
    "        H += g[qind]*(multiqop(oplist=oplist_g0)+multiqop(oplist=oplist_g1))\n",
    "        \n",
    "    #add bus\n",
    "    oplist_i = oplist_id.copy()\n",
    "    oplist_i[-1] = nop(buslevels) \n",
    "    H += fbus*multiqop(oplist=oplist_i)\n",
    "    return H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consider a single qubit coupled to the bus\n",
      "Shift in the qubit frequency from coupling to the bus: -3.491293 MHz\n",
      "Shift in the qubit frequency from equation: -3.500000 MHz\n",
      "\n",
      "\n",
      "Add Q1:\n",
      "Exchange Coupling: 3.564815 MHz\n",
      "Shift in the Q0 frequency from Q1 (on top of the bus shift): -0.233214 MHz\n",
      "Shift in the Q0 frequency from equation: -0.254819 MHz\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Hamiltonian parameters in MHz\n",
    "g = [70., 70.]\n",
    "fq = [5000., 5050]\n",
    "fbus = 6400.\n",
    "\n",
    "#Hamiltonian with 1 qubit coupled to the bus\n",
    "H1 = Hq(fq, fbus, [g[0], 0], qlevels=1, buslevels=1)\n",
    "\n",
    "print(\"Consider a single qubit coupled to the bus\")\n",
    "print(\"Shift in the qubit frequency from coupling to the bus: %f MHz\"%\n",
    "      (np.real(sorted(np.linalg.eig(H1)[0])[1])-fq[0]))\n",
    "print(\"Shift in the qubit frequency from equation: %f MHz\"%(qubit_dressed(fq[0], fbus, g[0])-fq[0]))\n",
    "print(\"\\n\")\n",
    "\n",
    "\n",
    "#Hamiltonian with 2 qubits coupled to the bus\n",
    "H1 = Hq(fq, fbus, g, qlevels=1, buslevels=1)\n",
    "Jcalc = J(fq[0], fq[1], fbus, g[0], g[1])\n",
    "print(\"Add Q1:\")\n",
    "print(\"Exchange Coupling: %f MHz\"%(np.abs(Jcalc)))\n",
    "print(\"Shift in the Q0 frequency from Q1 (on top of the bus shift): %f MHz\"%\n",
    "      (np.real(sorted(np.linalg.eig(H1)[0])[1])-qubit_dressed(fq[0], fbus, g[0])))\n",
    "print(\"Shift in the Q0 frequency from equation: %f MHz\"%\n",
    "      (qubit_dressed(qubit_dressed(fq[0], fbus, g[0]), qubit_dressed(fq[1], fbus, g[1]), Jcalc)-\n",
    "                                                        qubit_dressed(fq[0], fbus, g[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sect2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ZZ\n",
    "\n",
    "One of the consequences of the coupling is that there is a small ``ZZ`` term left over after the diagonalization. So \n",
    "\n",
    "$$H = \\sum_{i=0}^{N_{\\text{qubits}}} \\lambda_i \\hat{\\tilde{n}}_{i} $$\n",
    "\n",
    "is more realistically\n",
    "\n",
    "$$H = \\sum_{i=0}^{N_{\\text{qubits}}} \\lambda_i \\hat{\\tilde{n}}_{i} + \\sum_{i,j<i} \\xi_{ij} |ij\\rangle \\langle ij|$$\n",
    "\n",
    "where \n",
    "\n",
    "$$|ij\\rangle = \\hat{a}_i^{\\dagger} \\hat{a}_j^{\\dagger} |0\\rangle^N $$\n",
    "\n",
    "The ``ZZ`` is typically written in this way because the frequency of each qubit is empirically measured when the other qubits are in the ground state and so\n",
    "\n",
    "$$\\xi_{ij} = \\omega_{|ij\\rangle} - \\omega_i - \\omega_j $$\n",
    "\n",
    "Below we numerically compute ZZ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZZ: 137.097636 kHz\n"
     ]
    }
   ],
   "source": [
    "# Hamiltonian parameters in MHz\n",
    "g = [70., 75.]\n",
    "alpha = [-330, -330]\n",
    "fq = [5000., 5050.]\n",
    "fbus = 6400.\n",
    "\n",
    "#Hamiltonian with 1 qubit coupled to the bus\n",
    "H1 = Hq(fq, fbus, g, alpha=alpha, qlevels=2, buslevels=3)\n",
    "\n",
    "eigvals = np.real(np.linalg.eig(H1)[0])\n",
    "\n",
    "# need to be careful here because \n",
    "# of how to track from the bare states\n",
    "# this code works ok except near resonances\n",
    "inds = [0,0,0]\n",
    "bare_e = [fq[0], fq[1], fq[0]+fq[1]]\n",
    "closeness = [1000,1000,1000]\n",
    "for eind, eigval in enumerate(eigvals):\n",
    "    for i in range(3):\n",
    "        if np.abs(eigval-bare_e[i])<closeness[i]:\n",
    "            closeness[i] = np.abs(eigval-bare_e[i])\n",
    "            inds[i] = eind\n",
    "                       \n",
    "xi = eigvals[inds[2]]-eigvals[inds[1]]-eigvals[inds[0]]\n",
    "\n",
    "print(\"ZZ: %f kHz\"%(xi*1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sect3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dispersive Shifts\n",
    "\n",
    "When the qubit is coupled to a bus there is also a dispersive shift of the bus frequency depending on the state of the qubit. For coupling buses this is neglected (there should be no photons in the bus), however, this concept is the basis for qubit measurement with a separate set of readout buses. \n",
    "\n",
    "$$H = \\omega \\hat{n} + \\alpha \\frac{\\hat{n}}{2} (\\hat{n}-1)  + \\omega_{\\text{bus}} \\hat{n}_{\\text{bus}}+ g \\left(\\hat{a} \\hat{b}^{\\dagger} + \\hat{a}^{\\dagger} \\hat{b}\\right) $$\n",
    "\n",
    "is \n",
    "\n",
    "$$H = \\tilde{\\omega} \\hat{n} + \\tilde{\\omega}_{\\text{bus}} \\hat{n}_{\\text{bus}}+ 2\\chi \\hat{n} \\hat{n}_{\\text{bus}} $$\n",
    "\n",
    "where \n",
    "\n",
    "$$\\chi = \\frac{g^2}{\\omega-\\omega_{\\text{bus}}} \\frac{\\alpha}{\\omega+\\alpha-\\omega_{\\text{bus}}} $$\n",
    "\n",
    "We can look at this numerically as the difference of the cavity frequencies with the qubit in $|0\\rangle$ or $|1\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi: -238.359594 kHz\n",
      "Chi (Formula): -239.061142 kHz\n"
     ]
    }
   ],
   "source": [
    "# Hamiltonian parameters in MHz\n",
    "g = [50.]\n",
    "alpha = [-330]\n",
    "fq = [5000.]\n",
    "fbus = 6700.\n",
    "\n",
    "#Hamiltonian with 1 qubit coupled to the bus\n",
    "H1 = Hq(fq, fbus, g, alpha=alpha, qlevels=2, buslevels=3)\n",
    "\n",
    "eigvals = np.real(np.linalg.eig(H1)[0])\n",
    "\n",
    "# need to be careful here because \n",
    "# of how to track from the bare states\n",
    "# this code works ok except near resonances\n",
    "inds = [0,0,0]\n",
    "bare_e = [fq[0], fbus, fq[0]+fbus]\n",
    "closeness = [1000,1000,1000]\n",
    "for eind, eigval in enumerate(eigvals):\n",
    "    for i in range(3):\n",
    "        if np.abs(eigval-bare_e[i])<closeness[i]:\n",
    "            closeness[i] = np.abs(eigval-bare_e[i])\n",
    "            inds[i] = eind\n",
    "                       \n",
    "chi = (eigvals[inds[2]]-eigvals[inds[0]]-eigvals[inds[1]])/2\n",
    "\n",
    "print(\"Chi: %f kHz\"%(chi*1000))\n",
    "print(\"Chi (Formula): %f kHz\"%(Chi(fq[0],fbus,g[0],alpha[0])*1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Tags",
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
